/*

This file makes Figure 2

*/
global writeout 1 // change to 1 for writing out, 0 to not
global quietly quietly // comment out for noisy


$quietly { // top top quietly
  clear
  set more off
  set seed 1337
  set sortseed 1337

  capture cd "~/Desktop/RESTAT_CODE"

  global dir_root = "."
  global dir_code = "${dir_root}/code"
  global dir_data = "${dir_root}/data"
  global dir_output "${dir_root}/output"

  * Set graph font
    graph set window fontface "Times New Roman"

  // Written by table 3 code
  use "${dir_data}/horseracedata", clear

  sort state_fips date
  xtset state_fips date

  *scale variables for graphing
    gen pI_100k= pI * 100000
    gen prev_public2_100k = prev_public2_*100000

  *label variable
    label var pI_100k "Latent prevalence (Per 100,000)"
    label var prev_public2_100k "Reported prevalence (Per 100,000)"

  *new x-axis date2
    scalar sdate2 = td(24july2020)

  * Panel A
  preserve
    cap gen pI_pct = pI * 100
    gen prev_public2_pct = prev_public2_ * 100
    scalar limit_val = 1.2
    keep if date == td(01jul2020)

    gen jitter_bug = 0
    replace jitter_bug = -.02 if state2 == "TX"
    gen pI_pct2 = `=limit_val' + jitter_bug if pI_pct > `=limit_val' & !missing(pI_pct)

    gen stateup = state2+"^"

    graph twoway ///
        (scatter pI_pct prev_public2_pct if (pI_pct < `=limit_val'), mlabel(state2) msymbol(none)) ///
        (scatter pI_pct2 prev_public2_pct if (pI_pct >= `=limit_val'), mlabel(stateup) msymbol(none)) ///
        (lfit pI_pct pI_pct if (pI_pct <= `=limit_val'), leg(off)) ///
      , xtitle("Reported Prevalence in % (confirmed cases - deaths - recoveries)") ///
        ytitle("Latent Prevalence in %") ///
        ylabel(0 .5 1, nogrid) ///
        xlabel(0 .5 1) ///
        ysize(4) xsize(5.5) scale(1.2) scheme(s2color) graphregion(color(white))  ///
        name(state_scatter, replace)

    if $writeout graph export "${dir_output}\state_scatter_0703_smooth_color.pdf", ///
        replace name(state_scatter)
  restore


  ** Panel B

  *define dates for x-axis
    scalar sdate1 = td(1may2020)
    scalar sdate2 = td(01june2020)
    scalar sdate3 = td(24june2020)
    scalar sdate4 = td(24july2020)
    scalar sdate5 = td(1september2020)

    line prev_public2_100k date if state2=="UT" & date > mdy(04,15,2020) & date <= mdy(07,02,2021) ///
      , lwidth(thick) lcolor(navy) lpattern(solid) ///
    || ///
    line pI_100k date if state2=="UT" & date > mdy(04,15,2020) & date <= mdy(07,02,2021) ///
      , lwidth(thick) lcolor(cranberry) lpattern(dash) yaxis(2) ///
      xline(`=sdate3', lc(cranberry) lw(medthick) ) ///
      xline(`=sdate4',  lw(medthick) lc(navy)) ///
      xscale(range(22020 22098)) ///
      yscale(line)  ///
      xlabel(`=sdate1' `=sdate2' `=sdate3' `=sdate4' `=sdate5' , format(%tdMon-DD)) ///
      xtitle("") ytitle() ///
      legend(order(1 2) lab(2 "Latent prevalence") lab(1 "Reported Prevalence") region(color(white)) position(6))  ///
      ysize(4) xsize(5.5) scale(1.2) scheme(s2color) graphregion(color(white))  ///
      name(utah_line, replace)

  if $writeout graph export "${dir_output}/Utah_timeline_reported.pdf", replace name(utah_line)

} // end top top $quietlys
